!dk initsub
!*************************************************************
!*                                                           *
      subroutine initiv
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use cplot_com_M
      use boundary
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!*                                                           *
!*************************************************************
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer , dimension(8) :: ltab
      integer :: is, kkk, ireg
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double), dimension(8) :: dtab
      real(double), dimension(5) :: rint, eigint
      real(double), dimension(102) :: rmod, eigmod
      real(double) :: ri1, ri2, ri3, ri4, ri5, dcttp, del6,  angvsv, &
         drbar, gf, gbk, gl, gr, gt, gb, urght, ulft, vfrnt, vbck, wbtm, wtp, &
         omgf, omgbk, omgl, omgr, omgt, omgb, arg, rad, bleq23
      logical :: expnd, nomore, cntrun
!-----------------------------------------------
!
      namelist /generale/restart, ncyc_restart
      namelist /celin/precon, eandm, itmax, itsub, ibar, jbar, kbar, kplas, &
         explicit, scalsq, del1, del2, del3, del4, del5, del7, btorus, bvertcl&
         , efld0, q0, nphist, numgrid, taug, delta, iout, tmax, tmax0, tmax1, &
         tmax2, tmax3, testpar, lpr, ncyc2, pars, dcttd, dpttd, twfin, dt, &
         dcttd, dpttd, twfin, dto, cntr, dtoc, lama, lamb, lamc, adaptg, &
         newstuf, fields, rhoi, siei, beta, coll, collide, mgeom, rpl, rvac, &
         rwall, rmaj, dx, dy, dz, cartesian, plots, nfreq, nu_len, nu_comm, &
         smoothing, elle, nsmooth, el_rec, bx_ini, by_ini, bz_ini, eps_rec, &
         rnu_rec, epsil, lsplit, lcoal, npcel_control, mtresh, vtresh, &
         periodicx, periodicy, periodicz, wkl&
         , wkr, wkb, wkt, wke, wkf,ex_ext,ey_ext,ez_ext,forcefree, shock, &
         bcMx,bcMy,bcMz,avg,fluid_load,initial_equilibrium, ex_ini, ey_ini, ez_ini, &
         nome_file_input, ex_const, ey_const, ez_const, wave_amplitude, external_current, &
         relativity
!
      namelist /reg/nrg, npcelx, npcely, npcelz, xvi, yvi, zvi, drift, t_wall, &
         wsr_itg, el_itg, rhr, siepx, siepy, siepz, icoi, uvi, vvi, wvi, qom, &
         bcpl, bcpr, bcpb, bcpt, bcpe, bcpf, uvip, vvip, wvip, mpert, pert_gem, pert_x&
         , shear,inject_end,inject_front,inject_left,inject_right
!
      namelist /objects/nrg_obj, npcelx_obj, npcely_obj, npcelz_obj, chii_obj, &
         rho_obj, icoi_obj, theta1_obj, theta2_obj, phi1_obj, phi2_obj, &
         r_minor1_obj, r_minor2_obj
!l    *** additional declarations from initiv ***
!
      equivalence (rint(1), ri1), (rint(2), ri2), (rint(3), ri3), (rint(4), ri4&
         ), (rint(5), ri5)
      write (*, *) ' initiv: reading header'
      read (5, 1111) nome_run
 1111 format(a80)
      read (5, nml=generale)
      write (59, nml=generale)
      if (restart) return
!      if (restart) go to 9527
!
      xvi=0d0
      yvi=0d0
      zvi=0d0
      do ireg=1,nreg
      xvi(1:8,ireg)=(/0.,1.e10,1.e10,0.,0.,1.e10,1.e10,0./)
      yvi(1:8,ireg)=(/0.,0.,1.e10,1.e10,0.,0.,1.e10,1.e10/)
      enddo
      pert_x=0.d0
      pert_gem=0.d0
      shock=.false.
      forcefree=.false.
      ex_ext=0d0
      ey_ext=0d0
      ez_ext=0d0
      ex_ini=0d0
      ey_ini=0d0
      ez_ini=0d0
      ex_const=0d0
      ey_const=0d0
      ez_const=0d0
      cntrun=.false.
      iout = 0
      dto = 0.
      dtoc = 0.
      i = 11
      touch = .FALSE.
      dcttd = 0.
      dcttp = 0.
      dpttd = 0.
      tmax0 = 0.
      tmax1 = 0.
      tmax2 = 0.
      tmax3 = 0.
      kplas = 1
      del1 = 0.
      del3 = 0.
      del4 = 0.
      del5 = 0.
      del6 = 0.
      del7 = 0.
      nay = lay
      nay2 = 1
      mgeom = 1
      nrg = 2
      drift(:nrg) = .FALSE.
      testpar = .FALSE.
      periodic = .TRUE.
      collide = .FALSE.
      lama = 0.
      lamb = 0.
      lamc = 0.
      adaptg = .FALSE.
      fields = .TRUE.
      plost_neg = 0.
      plost_pos = 0.
      bcpb = 2
      bcpt = 2
      external_current = .FALSE.
      wave_amplitude = 0d0
      relativity = .FALSE.
      periodicx = .TRUE.
      periodicy = .TRUE.
      periodicz = .TRUE.
      bc_vtoctov = .false.
!
!	Parametres not changed from case to case
!
!
      divergence_cleaning=.true.
      pars=.5
      nphist=001
      cartesian=.true.
      smoothing=.false.
      cntr=1.
      plots=.true.
      nfreq=1
 	  del2=.001
      del3=.00
      del4=.000
      del5=+.000
      mgeom=11
      lpr=1
      SCALSQ=0.0
      taug=0.20
      numgrid=50
      explicit=.false.
      eandm=.true.
      adaptg=.false.
      newstuf=.true.
      fields=.true.
      delta=0.95
      dtoc=100000.
      beta=1.
      siei=.00
      rhoi=1.
      btorus=0.
      q0=0.00
      efld0=0.
      bvertcl=0.00
      tmax=3.
      rpl=1.
      rvac=1.
      rwall=1.
      rmaj=0.0
!
!     default plot interval
!
      t = 0.0
      twfin = 0.
      ixto = 1
      if (dto(ixto) == 0.) dto(ixto) = twfin/4.0
      tout = t + dto(ixto)
      call date_and_time(ddate,hour)
      write(*,*)'date (ccyymmdd) = ',ddate
      write(*,*)'time (hhmmss.sss) = ',hour
9527  continue
!
!ll   1.  read parameters from cards:
!         --------------------------
      write (0, *) 'initiv: reading celin'
      read (5, nml=celin)
      write (59, nml=celin)

      if(ncyc2*dt/dto(1)>150.d0) then
	write(*,*) 'too many output frames, increase dto'
	stop
      endif
!l    (rmod is a factor of 10 too large on input file -- adjust)
      rmod(:13) = rmod(:13)/10.
      i = 14
!
!
!
!ll   2.  check if continuation of previous run:
!         -------------------------------------
      if (ibar<0 .and. jbar==0) stop
      if (ibar == 0) then
!
!l    2.1 continuation -- read parameters and all data from tape:
!         ------------------------------------------------------
         cntrun = .TRUE.
!
         omg = angv
         if (t < tangv) omg = t*omgdt
!
         angvsv = angv
!
!
         if (angv /= angvsv) then
            tangv = amax1(tangv,t + dt)
            omgdt = (angv - omg)/(t - tangv)
         endif
!
         ncyc1 = ncyc + 1
!
      else
!
!l    2.2 read further parameters from data cards:
!         ---------------------------------------
         write (0, *) 'initiv: reading reg'
         read (5, nml=reg)
         if (nrg > nsp) then
            write (59, 9901)
 9901       format(" nrg is larger than nsp")
!      stop
         endif
!
         anydrift = .FALSE.
         do is = 1, nrg
            if (.not.drift(is)) cycle
            anydrift = .TRUE.
         end do
!
         strait = 0.0
         toroid = 1.0
         if (rmaj >= 1.0E+05) rmaj = 0.
         if (rmaj == 0.0) then
            strait = 1.0
            toroid = 0.0
         endif
!
         kvac = kplas + 2
         kplas = kbar
         if (.not.cartesian) then
            dx = rwall
            dy = rmaj
         endif
         wsa = rpl
         a = rvac
!
         ncyc1 = 1
         if (ncyc2 < 0) ncyc2 = 77777
      endif
!
      write (*, *) 'reading objects'
      read (5, nml=objects)
!
!
      if (tmax0 <= 0.) tmax0 = 1.
      if (tmax1 <= 0.) tmax1 = 1.
      if (tmax2 <= 0.) tmax2 = 1.
      if (tmax3 <= 0.) tmax3 = 1.
      if (tmax <= 0.) tmax = (tmax0 + tmax1 + tmax2)/3.
!
!
!ll   3.  print input data
!         ----------------
      kpr = kpt
!
      write (59, 1111) nome_run
      write (59, nml=celin)
      write (59, nml=reg)
      kpr = kfm
!
      if (.not.cntrun) then
!
!
!ll   4.  calculate data for run-setup:
!         -----------------------------
         ibp1 = ibar + 1
         jbp1 = jbar + 1
         kbp1 = kbar + 1
         ibp2 = ibar + 2
         jbp2 = jbar + 2
         kbp2 = kbar + 2
         kvac = kplas + 2
         kvp = kvac + 1
         kvm = kvac - 1
         kvm2 = kvac - 2
         iwid = lay
         jwid = ibp2*iwid
         kwid = jbp2*jwid
         iper = ibar*iwid
         jper = jbar*jwid
         fibar = float(ibar)
         rfibar = 1./fibar
         fjbar = float(jbar)
         rfjbar = 1./fjbar
         fkbar = float(kbar)
         rfkbar = 1./fkbar
         rbjbkb = 1./(ibar*jbar*kbar)
         pi = acos(-1.d0)
         t = 0.
         mu(1) = 0.
         lam(1) = 0.
         gm1 = 0.6666667
         ixto = 1
         tout = t + dto(1)
         pttd = t + dpttd
         ncyc = 0
         numtd = 0
         drbar = 0.9*rwall
!      dtpos=dt*0.1
!      dt=dtpos
         dtc = dt
         dtv = dt
         nomnit = 10
         numit = 10
         rdt = 1./(dt + 1.E-20)
         anc = 0.
         anc = 0.05
         kxi = -1
         colamu = (1. + 1.67)/(lam(1)+mu(1)+mu(1)+1.E-10)
         omanc = 1. - anc
         lam(2) = lam(1)
         mu(2) = mu(1)
         gf = 0.
         gbk = 0.
         gl = 0.
         gr = 0.
         gt = 0.
         gb = 0.
         gf = 1.
         gbk = 1.
         gl = 1.
         gr = 1.
         gt = 1.
         gb = 1.
         urght = 0.
         ulft = 0.
         vfrnt = 1.
         vbck = 1.
         vfrnt = 0.
         vbck = 0.
         wbtm = 0.
         wtp = 0.
         omgf = 1. - gf
         omgbk = 1. - gbk
         omgl = 1. - gl
         omgr = 1. - gr
         omgt = 1. - gt
         omgb = 1. - gb
         thrd = 1.0/3.0
         sxth = 1.0/6.0
         itrlmp = 25
!
!     parameters for adaptive regriding
!
         frrg = 0.01
         alrg = 1.
         dtrg = dt
         trgrd = t
         trglst = t - dt
!
!     define periodic index switch
!
         irt(1) = 1
         irt(2) = 1
         irt(3) = 0
         irt(4) = 0
         irt(5) = 1
         irt(6) = 1
         irt(7) = 0
         irt(8) = 0
         ibk(1) = 0
         ibk(2) = 1
         ibk(3) = 1
         ibk(4) = 0
         ibk(5) = 0
         ibk(6) = 1
         ibk(7) = 1
         ibk(8) = 0
!
         h = 2.0*pi/dz
         wghtj = btorus*dz*rfjbar
!
         tflux = 0.5*fibar*sin(2.*pi*rfibar)*(rwall*rwall - rvac*rvac)
!
!     initial poloidal inductance estimated assuming coaxial cylinders
         elpp = 0.5*dz*log(rwall/rvac)/pi
      endif
!
      wsa = amin1(rpl,rvac)
!
      if (q0 == 0.) then
         rq0 = 0.
!
!        parameters for scyllac field
!
         arg = h*rwall
         rad = 0.01
!
         do kkk = 1, 4
            if (kkk == 2) rad = wsa
            if (kkk == 3) rad = rvac
            if (kkk == 4) rad = rwall
            arg = h*rad
!            del0= rad*dl0
         end do
!
      else
!
         rq0 = 1./q0
         bl0 = 0.
         bl0sv = 0.
         bl1 = 0.
         bl1sv = 0.
         bl2 = 0.
         bl2sv = 0.
         bl3 = 0.
         bl3sv = 0.
         bvrt = 0.
         bvrtsv = 0.
         bleq23 = 0.
         bl23 = 0.
         rdl0sv = 0.
         rdl1sv = 0.
         del2sv = 0.
      endif
!
!
!ll   5.  define mesh and initialize diagnostics:
!         ---------------------------------------
      call celindex (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, kwid, ijkcell, &
         ijkctmp, ncells, ijkvtx, nvtxkm, nvtx)
      write(*,*)'initiv,celindex',ibp1,jbp1,kbp1,ncells,cntrun
!
      if (.not.cntrun) then
!
!l    5.1 calculate mesh:
!         --------------
         call mshset
      endif
!
!
      return
      end subroutine initiv
